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I.  Introduction 


Modeling  of  nonequilibirum  processes  in  a  low-temperature  partially  ionized  plasma  is  of 
particular  interest  to  a  wide  range  of  technical  fields  such  as  gas  discharge,  electric  propulsion, 
spectroscopic  and  laser  diagnostics,  and  material  science1-3.  The  complexity  of  the  model  is 
largely  due  to  the  characterization  of  various  collisional  and  radiative  processes,  occurring  at 
a  wide  range  of  spatial  and  temporal  scales4,5.  Although  the  fundamental  physical  processes 
may  be  individually  known,  it  is  not  always  clear  how  their  combination  affects  the  overall 
operation,  or  at  what  level  of  detail  this  process  needs  to  be  modeled.  The  current  state  of 
the  art  for  modeling  detailed  chemical  kinetics  of  a  low  temperature  plasma  is  the  collisional- 
radiative  (CR)  model,  first  proposed  by  Bates  et  al.  in  19626,7.  CR  models  are  now  commonly 
used  in  studies  of  plasma  discharge,  plasma-assisted  combustion,  and  hypersonics8-13.  The 
advantage  of  a  CR  model  is  two- fold.  First,  strong  deviations  from  equilibrium  of  the  internal 
states  can  be  captured  accurately  when  CR  models  are  employed.  In  addition,  ab  initio 
cross  section  data  can  be  directly  incorporated  in  the  CR  model,  leading  to  a  very  accurate 
prediction  of  the  thermochemical  kinetics  of  the  system. 

There  are,  however,  two  issues  with  this  modeling  approach.  The  first  arises  from  the 
complexity  of  the  physical  processes  needed  to  be  captured  in  the  model.  The  required 
level  of  detail  of  the  CR  model  is  typically  not  known  a  priori  and  is  possibly  changing  in 
a  dynamical  fashion  as  the  system  evolves  in  time.  This  can  be  resolved  by  coarse-graining 
techniques,  which  reduce  the  complexity  of  the  kinetics  to  avoid  solving  a  large  system  of 
equations14-16.  The  second  issue  comes  from  translational  nonequilibrium,  often  found  in  a 
discharge,  where  we  have  both  a  bulk  plasma  (continuum)  and  a  highly  energetic  component 
(kinetic),  e.g.,  electrons  emitted  from  the  cathode.  A  proper  treatment  of  this  energetic 
beam-like  component  requires  extending  the  solution  of  the  CR  kinetics  to  the  so-called 
non-Maxwelllian  regime17-19.  These  simulations  are  typically  very  expensive  and  therefore 
limited  to  zero-  or  one- dimensional  systems.  In  addition,  space  charge  effects  within  the  bulk 
plasma  can  also  become  important,  requiring  further  separation,  i.e.,  ions  and  electrons.  A 
natural  solution  to  this  problem  is  to  use  a  hybrid  approach,  decomposing  the  system  into 
a  continuum  and  a  kinetic  components20.  The  continuum  component  can  be  solved  by  fluid 
equations,  and  the  kinetic  component  by  (for  example)  a  particle  method.  Although  the 
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idea  seems  quite  intuitive,  proper  treatment  of  the  coupling  between  the  fluid  and  kinetic 
components  is  highly  non-trivial.  There  appears  to  be  no  unique  and  consistent  coupling 
methodology,  and  the  choice  is  highly  problem-dependent.  This  coupling  issue  is  current 
being  addressed  by  many  researchers  and  is  outside  of  the  scope  of  the  current  work. 

We  are  concerned  here  with  an  alternative  approach,  the  so-called  multifluid  model,  which 
decomposes  the  plasma  into  several  fluid  components.  For  example,  in  a  discharge  configu¬ 
ration,  one  could  have  4  different  fluids,  namely  the  neutrals,  ions,  bulk  electrons,  and  the 
energetic  electrons.  The  only  required  assumption  is  that  collisions  among  particles  within 
the  same  fluid  are  sufficiently  fast  to  maintain  a  Maxwellian  distribution.  The  validity  of 
this  assumption  is  not  always  well  known.  Nevertheless  this  approach  is  attractive,  since 
it  is  much  faster  than  a  fully  kinetic  solver,  and  unambiguous  since  at  the  same  time,  the 
approach  relies  on  kinetic  theory  for  the  treatment  of  coupling  terms  between  different  fluids. 
Furthermore,  one  can  rely  on  fast  implicit  methods  to  solve  these  fluid  equations  and  exam¬ 
ine  long  time  behavior  of  the  system;  this  offers  a  big  advantage  over  kinetic  simulations, 
which  are  only  suitable  for  problems  with  short  time  scales.  Multifluid  models  are  commonly 
used  in  simulations  of  astrophysical  plasmas  (see  for  example21,22). 

The  most  classical  work  in  multifluid  plasma  modeling  is  due  to  Braginskii23,  who  derived 
fluid  equations  for  a  fully  ionized  plasma,  using  a  Chapman- Enskog  closure.  Braginskii’s  work 
has  been  successively  refined  by  several  authors,  with  particular  emphases  on  improving  the 
transport  coefficients  and/or  including  interaction  with  neutral  species24-2'.  Burgers,  on  the 
other  hand,  presents  a  rather  general  framework  for  the  modeling  of  elastic  collisions28.  These 
include  both  neutral  collisions  and  charged  particle  collisions;  the  methodology  is  applicable 
for  a  general  system  of  moment  equations  beyond  the  standard  five-moment  model.  Burgers 
also  introduces  a  simplified  model  for  reactive  collisions  using  a  Bhatnagar-Gross-Krook 
(BGK)  collision  operator. 

In  the  current  work,  we  present  a  self-consistent  model  for  inelastic  collisions  within  the 
multifluid  framework.  The  model  is  derived  from  kinetic  theory  and  obeys  the  principle  of 
detailed  balance  (DB),  which  we  show  to  be  an  essential  property  to  ensure  that  the  system 
approaches  the  correct  equilibrium  limit.  We  focus  on  characterizing  the  exchange  source 
terms  due  to  collision,  namely  mass,  momentum  and  energy  exchanges  (the  hydrodynamics 
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and  transport  fluxes  can  be  added  following23  when  considering  a  non-uniform  plasma).  We 
will  show  that  in  most  cases  none  of  these  terms  can  be  neglected,  and  they  have  complex 
dependencies  to  microscopic  quantities  of  the  interaction,  e.g.,  multiply-differentiated  cross 
sections.  We  present  the  general  description  of  the  collision  kinematics  and  derive  the  ex¬ 
change  terms  for  the  case  of  excitation  and  deexcitation  processes.  Although  we  are  mostly 
interested  here  in  electron- impact  collisions  and  atomic  transitions,  we  keep  the  formulation 
as  general  as  possible,  such  that  the  application  to  other  species  and  chemistry  (e.g.  proton- 
impact,  molecular  vibrational  transitions,  charge-exchange,  etc.)  is  a  straight-forward  exten¬ 
sion.  The  case  of  ionization  and  recombination,  and  other  three-body  processes,  is  currently 
under  examination,  using  the  same  basic  formulation  presented  here. 

The  rest  of  the  paper  is  organized  as  follows.  The  derivation  of  the  exchange  source  terms 
is  given  in  Sec.  II,  by  first  introducing  the  description  of  the  transfer  integral,  and  then 
present  ingthe  derivation  of  the  exchange  rates  in  the  following  subsections.  In  Sec.  Ill, 
we  show  the  numerical  evaluation  of  the  multifluid  rates,  verify  the  results  with  Monte 
Carlo  calculations,  and  perform  zero-dimensional  calculations  utilizing  the  multifluid  rates. 
Finally,  conclusions  and  a  summary  of  the  present  findings  are  given  in  Sec.  IV.  We  also 
provide  several  appendices  to  elaborate  on  the  derivation  of  the  exchange  source  terms  and 
the  description  of  the  numerical  simulation. 

II.  Rate  Derivation 
A.  Transfer  integral 

Let  us  consider  an  inelastic  collision  between  two  particles  s  and  t,  such  that  the  particle  t 
changes  its  internal  state.  The  particles  s  and  t  are  respectively  the  scattered  and  target  in 
the  laboratory  frame  of  reference  (LAB).  The  former  will  be  identified  as  the  electron  and 
the  target  as  the  atom,  but  we  will  keep  the  general  s,  t  notation  until  explicit  assumptions 
and  approximations  are  made,  such  as  neglecting  terms  of  the  order  of  the  mass  ratio  ms/mt 
for  final  expressions.  Following  Appendix  A,  the  initial  velocities  are  vs,  vt,  where  v  =  u+c 
and  u  is  the  fluid  mean  velocity  in  the  LAB  frame,  and  post-collision  values  are  indicated 
by  a  prime,  i.e. : 


s(vs)  +  t(vt)  ->•  s'(v's)  +  t'(v'i) 


(1) 
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We  make  here  two  assumptions:  1)  the  collision  produces  only  two  particles,  which  may  or 
may  not  belong  to  the  same  fluids  as  the  initial  reactants,  and  2)  the  masses  of  individual 
particles  are  the  same  before  and  after  the  collision,  e.g.  m's  =  ms,  such  that  mass  conserva¬ 
tion  is  automatically  obtained.  Both  of  these  will  be  revisited  in  a  follow-on  paper,  dealing 
with  ionization  and  recombination.  Defining  the  energy  transfer  to  and  from  the  internal 
modes  to  be  represented  by  Ae,  we  have  the  following  energy  conservation  constraint  on  the 
relative  velocity  g  where  g  =  vs  —  vf  (see  Appendix  A): 


g2  =  g'2  + 


2Ae 

h 


(2) 


For  excitation,  the  transferred  energy  is  a  positive  and  fixed  value  Ae  =  £*,  the  energy 
gap  between  the  levels,  while  for  ionization  it  is  a  continuum  of  values:  Ae  E  [£*,£],  where 
£  =  |/Ug2  is  the  available  kinetic  energy  in  the  center-of-mass  (COM)  frame.  In  the  limit 
A£  — >  0,  the  collision  is  elastic.  We  will  keep  the  same  relations  for  the  reverse  process,  for 
which  the  primed  variables  are  post-collision  and  non-primed  refer  to  pre-collision,  such  that 
for  deexcitation,  Ae=—e*. 


We  can  then  define  a  transfer  integral  of  the  collision  operator  between  the  two  species  s 
and  t28. 


qfst  =  nsnt  /  d3vsd3vt  fsft  g  /  0  du(vs,  vt;  v's,  v't) 


(3) 


where  g  is  the  magnitude  of  the  relative  velocity  (g  =  |g|),  du  is  the  differential  cross  sec¬ 
tion  (DCS),  and  0  is  any  moment  variable  exchanged  during  the  collision.  We  now  follow 
Appendix  B,  starting  with  the  following  transformations: 


V* 

g 

and  7 


V  -  U  +  7g 

g  -  w 

riTt-TB) 

msTt+mtTs 


,*  =  MTsTt  2  =  2kT* 

rnsTt  +rntTs  °  M 

~  rnsTt+rntTs  2  2  kT 

M  u 


(4a) 

(4b) 

(4c) 


where  the  relative  mean  velocity  w  =  us  —  ut.  The  product  of  the  two  Maxwellian  distribu¬ 
tions  fs  ■  ft  is  expressed  in  terms  of  the  product  of  two  other  Maxwellians,  fy*  ■  fg,  for  the 
COM  velocity  and  relative  velocity  respectively.  Inserting  (B.14-B.16)  into  (3),  the  transfer 
integral  can  be  written  as  follows: 


dOi  =  nsnt- 


—  s 

2  /7V 


d3\* 


V*2/a2 


7T2a 


—  s 

9  A/O 


d3 ge  g2/a2g  /  fjduj( g;  g') 


irza 


(5) 
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Figure  1.  Frame  rotation  and  relative  orientation  of  (a)  w  and  g  and  (b)  g  and  g'.  The  rotation 
operator  matrix  R(<p,  9)  (or  R(p ,  y))  is  defined  such  that  g  =  R(ip,  9)  ■  w  and  g'  =  R(p ,  y)  •  g- 


Note  that  in  the  COM  reference  frame,  the  DCS  only  depends  on  the  relative  velocities,  i.e., 
duj(vs,  vt]  Vs,  v't)  =  du(g]  g'),  and  can  be  expressed  as: 

duj(g]g,')  =  <jst(g,n')dn'  (6) 


where  O'  is  the  solid  angle  between  the  initial  and  final  relative  velocities,  i.e.,  diV  —  dp dcos y 
with  g-g'  =  gg1  cos  y.  Without  loss  of  generality,  we  can  now  choose  a  reference  frame  (LAB) 
such  that  the  relative  mean  velocity  w  is  aligned  with  the  z  axis,  as  shown  in  Figure  1.  Thus, 
the  unit  vectors  g,  g'  are  obtained  by  subsequent  rotations  of  the  (x,  y,  z)  frame.  Using  the 
abbreviated  notation  cv  =  cos  99,  s ^  =  sin  tp,  etc,  we  dehne  this  rotation  operator  by  the 
matrix: 


c^ce  S(p  CtpSo  ^ 


R(<P,0)  = 


S(pC@  C(p  S^So 

-so  0  Co 


and  g  =  —  =  R(<p,  0)  ■  z  = 
9 


(  CvSo  ^ 
SipSQ 

Ce  ! 


(7) 


Similarly,  the  post-collision  relative  velocity  is  rotated  by  the  angles  (p,  x),  such  that  g'  = 

R(p,x)  ■  g- 


Using  d3g  =  g2dgdipdco1  and  equation  (6),  the  transfer  integral  can  be  written  as: 

^ st  =  ^±e~w2'a2  •  J  d3\*fv*  •  J  dgg3e~a2/°2- 

dipdco  e2s™C0/“2  f  dpdcxil)(jst{g,Xl') 


7T2CT 


(8) 
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Let  us  now  assume  that  the  moment  variable  can  be  expanded  in  terms  of  powers  of  V*: 


0  =  a  +  bV*  +  cV*2  +  ... 


where  a,b,c...  are  functions  of  the  remaining  velocity  variables,  and  let  us  perform  the 
integration  over  V*.  Note  that  we  have: 


d3V*fv*  =  1 


d3\*\*fy,  =  0 


the  latter  by  reasons  of  symmetry.  Thus,  as  long  as  0  does  not  contain  terms  quadratic  (or 
higher)  in  V*,  a  condition  satisfied  throughout  this  work,  we  can  eliminate  the  integration 
over  V*,  keeping  only  the  terms  which  are  independent  of  V*.  Also  by  symmetry,  the  DCS 
a st  does  not  depend  on  the  angle  p,  and  we  can  write: 


a.t(g,M)  =  <?,t(g)-Q(g,x)  S.t.  /  dp  dcxQ(g,  x)  =  1 


(9) 


where  Q  is  the  angular-dependent  DCS.  More  generally,  we  will  define  the  averaging  of  any 
function  0  over  the  scattering  angles  as: 


■+i 


W')„.  =  2>r  /  <fc*V>S(9,x) 


(10) 


'-i 


A  trivial  integration  over  ip  yields: 

=  / dgg3e-92'a2ast(g )•  [+'dcee2awc°'a2  (0) 


7T  2  Q?3 


Q' 


(ii) 


We  now  define  the  following,  normalized  energy  variables, 


z  = 


A  = 


1  2 
2 


z  =  max  0, 


Ae 


(12) 


kT  kT  V  kT, 

where  T  is  dehned  in  Appendix  B.  Using  g3dg  =  2ede/p2  and  a  further  change  of  variables, 
we  finally  obtain: 


8  kT  _A 

Vst  =  nsnt  I  -  e 


0  +  1 


izze~2Wlt(z)  -J  dce2^  • 


(13) 


with  gf  a  thermal  velocity  based  on  the  average  temperature  T.  Note  that  we  have  left  the 
variable  -0  undetermined,  and  since  it  could  potentially  depend  on  all  integration  variables 
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(z,x,8),  V’  must  be  kept  inside  all  integrals.  We  will  see  what  simplifications  can  be  made 
next,  depending  on  which  moment  variables  we  are  integrating. 

The  lower  limit  of  integration,  z*,  is  zero  for  elastic  collisions  or  for  exothermic  reactions. 
Thus,  equation  (13)  is  a  general  formula,  which  applies  equally  to  excitation  (z*  >  0)  and 
deexcitation  ( z *  =  0)38.  Let  us  first  consider  an  excitation  collision: 

s(vs)  +  t(E£,  -»•  s(v's)  +  t(Eu,  v't) 

where  t  and  u  denote  the  lower  and  upper  energy  states,  respectively.  We  made  here  the 
assumption  that  both  states  (£,  u)  belong  to  the  same  fluid,  so  that  the  particle  indices  (s,  t) 
are  kept  the  same,  but  this  is  not  necessary.  From  eq.  (2),  energy  conservation  implies 
Ae  =  Eu— Eg  >  0.  We  can  then  define  normalized  energy  variables  for  this  case,  x,x'  and  x*: 

8* 

x*  =  ~  >0;  x  =  x  —  x*  >  0  (14) 

kT 

Note  that  x  (x')  is  the  normalized  kinetic  energy  of  the  initial  (final)  products  of  excitation 
respectively,  and  that  x*  is  the  normalized  energy  threshold,  always  a  positive  quantity.  For 
excitation,  we  can  use  (13)  with  the  following  identifications: 

z  =  x,  z*  =  x*,  z  =  x' ,  and  nt  =  nt  (15) 

For  a  deexcitation  collision,  we  have  the  reverse  (u  — >  £),  i.e.  Ae  <  0: 

s(vs)  +  t(Eu,  vf)  ->•  s(v's)  +  t(Ee ,  v't)  (16) 

Equation  (13)  is  again  still  valid,  if  we  now  make  the  following  ident ideations: 

z  =  x\  z*  =  0,  z'  =  x,  and  nt  =  nu  (17) 

Therefore,  in  all  cases  the  variable  always  x  refers  to  the  larger  kinetic  energy  (before  ex¬ 
citation  or  after  deexcitation)  and  x'  refers  to  the  smaller  value  (after  excitation  or  before 
deexcitation).  We  can  therefore  define  two  cases  of  (13): 

roo  -i  p  +  1 

'k li  =  nsnegfe~x  dxxe~xale(x)--  dcee2^Ce  ■  {ip)n,  (18a) 

Jx*  2  J_1 

=  nsnugfe~x  dx' x' e~x'  ^u(x')  •  ^  J  dcee2V^Ce  ■  (^)n,  (18b) 


where  the  superscripts  f,  4-  indicate  excitation  and  deexcitation  respectively  (note  the  change 
of  subscript  from  st  to  si  for  excitation,  and  su  for  deexcitation).  It  is  worth  pointing  out  that 
the  averaging  over  the  scattering  angle,  i.e,  (0)q/,  has  to  be  done  with  the  corresponding 
angular-dependent  DCS  Q,  e.g.,  (V’)n'  =  27 t  f  cxyjQl(  for  excitation.  However,  from  time 
reversal  we  have  Q\g  —  G\ui  so  for  simplicity,  we  do  not  differentiate  (0)  between  excitation 
and  deexcitation.  Note  that  the  integration  over  the  6  angle  remains  to  be  performed.  If 
the  moment  variable  0  can  be  expanded  in  terms  of  power  of  cos  9,  we  can  then  define  the 
following  set  of  functions: 

C(fc)(0=A4  [+1dyyke2^  (19) 


i- i 


where  A4  is  a  normalizing  factor.  In  particular,  we  have: 
1 


c(0)«)  = 


r'dye^  =  Sillh(2?) 

/-i  25 

f  +  1  O 

,2 iy  _  6 


=  -&:  j^yy^  =  -^2 


cosh (24)  — 


sinh(2£) 

2e 


s.t.:  lim  C(0)  =  1 

(20a) 

s.t.:  lim  C(1)  =  1 
£->o 

(20b) 

In  the  CR  model,  each  internal  state  is  treated  as  a  pseudo-species,  so  the  rate  of  change  in 
number  density  for  each  state  (ng,  nu)  is  taken  into  account  separately.  We  can  now  examine 
the  specific  form  taken  by  the  transfer  integral,  according  to  the  chosen  moment  variable, 
starting  from  (8),  (13),  or  (18). 


B.  Zero*h-order  moment:  number  density 

The  rate  of  change  of  the  number  density  due  to  an  inelastic  collision  of  type  (1)  can  be 
obtained  by  setting  0  =  1  in  (13),  so  the  average  over  all  the  scattering  angle  is  trivially 
removed: 

/oo  ^  p+1 

dz  z  e~z  a st{z)  •-  J  dcge2y/TzC9  (21) 

The  integration  over  dcg  yields  the  function  defined  in  eq.  (20).  We  can  now  express  the 
rates  for  transitions  between  two  atomic  levels  i,  u,  by  making  the  appropriate  substitutions 
for  the  energy  variables.  For  the  case  of  an  excitation  [i  — y  u),  and  according  to  eq.  (15), 
we  define  the  variable  x  as  the  normalized  kinetic  energy  of  the  reactants  (s,  t)  in  the  COM 
frame,  prior  to  the  collision:  therefore  in  this  case,  z  =  x,  z*  =  x*  >  0  and  nt  =  ng.  Thus, 

roo 

=  nsnigfe~x  /  dx  x  e~xale(x)(('0\VXx)  (22) 

Jx* 
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For  deexcitation  (u  — >■  £),  the  rate  of  change  of  number  density  follows  from  (17): 


Tisu  =  nsnugfe  x  /  dx' x' e  x'  a\u{x')C{0\V\x') 


Both  of  these  quantities  are  positive,  hence  the  resultant  rates  equations  are: 


(23) 


due  _  .p-f- 

~dt~  s*~ 


dnu 

dt 


and 


dni 

dt 


=  +Tl  = 

1  ^  shi. 


dnu 

dt 


In  the  case  of  electron- impact  processes  (s  =  e),  we  can  neglect  terms  of  order  me/M,  and 
for  an  atomic  transition  between  levels  £—)•«,  we  obtain: 


frf  =  HeV'e  e  A  /  dxxe  X  CT^x)  (V\x) 


(24) 


where  ve  =  \J In  the  limit  of  thermal  plasma  where  multifluid  effects  are  weak,  i.e. 
A  — >  0,  we  obtain: 

POO 

F le  —  nen£Ue  /  dxxe~xa\^x)  (25) 

Jx* 

which  is  exactly  the  expected  result  for  a  single-fluid  plasma. 

Using  the  Klein-Rosseland  relation  for  detailed  balance30, 

=  °L(x')x'gu  (26) 

where  tfcyjju  are  the  degeneracies  of  the  lower  and  upper  atomic  levels  respectively,  we  can 
write  the  excitation  rate  as  follows: 


=  nsntgf  e  x—e  x* 


dx' x' e  x’  C^0)(A/A(x*+a;/))  o\u{x') 


One  can  then  easily  extract  reaction  rates,  for  example: 


rL  =  ™  si  ■  nsnl 


r t  ■ 

an  i  an  i 


It  is  instructive  to  consider  the  ratio  of  these  rates: 

[n  1  f™dx'  x'  e~x'  C (  a/A  (x1 + x* ) )  alu  {x> ) 


i 

su 


— e-a 


(27) 


(28a) 

(28b) 


(29) 


J^°dx'  x1  e~x'(^  (Vl\x')aju(x') 

The  first  term  in  brackets  is  the  traditional  Boltzmann  equilibrium  relation;  the  second 
term  contains  the  correction  due  to  the  multifluid  effects,  and  appears  only  through  the 
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function  (20-a).  A  Taylor  expansion  near  A  =  0  yields  (with  an  obvious  definition  of  the 
Boltzmann  function  B): 


w 


si 


l&su 


-x* 

L  m 

Beu(T) 


f^dx'x'e  x' 


1  + 


2A(a:*-Hr') 


tL(t') 


f™dx'x'e~x'  [1  +  ^f\  <r^su{x') 

,  2Ax*  ..2, 

H — - 1-  o(A2) 


(30) 


Thus,  we  recover  the  expression  for  Boltzmann  equilibrium  in  the  thermal  (single-fluid)  limit 
(A  — »  0).  Note  that  the  correction  term  increases  with  the  energy  threshold,  i.e.  transitions 
between  high  levels  (x*  — >  0)  will  not  be  affected  very  much  by  the  multifluid  effects,  while 
the  impact  will  be  stronger  for  excitation  from  low  energy  levels,  with  high  energy  gaps.  For 
elastic  collisions  (x*  =  0),  the  ratio  of  rates  is  exactly  given  by  the  ratio  of  degeneracies. 


C.  First-order  moment:  momentum  density 

Consider  now  the  forward  reaction  (1)  and  the  corresponding  loss  of  momentum  to  the 
particles  with  velocity  v,s.  The  transfer  variable  in  this  case  is  ^  =  msv„  and  starting  from 
equation  (8),  the  contribution  to  the  momentum  equation  is: 

R-7  =  -- •  [  d3V*fv *  •  [  dgg3  e~g2/a2ost{g)  ■  f  d(pdce  e29WCe/a2  (msvs)n,  (31) 
7T2CT  J  J  J 

Similarly,  the  gain  in  momentum  is  given  by  the  production  of  new  particles  with  velocity 
v's: 


R^  =  +-¥^:-  [  d3V*fv*-  f  dgg3e  9'/a2  crst(g)  ■  [  dtpdcee29WCe/a2  (msVs)n,  (32) 
7T2CT  J  J  J 

Using  the  relation: 

ms(vs  —  v's)  =  p(g  —  g')  (33) 

we  verify  that  the  integrand  does  not  depend  on  V*  and  its  integration  is  trivially  removed. 
The  net  rate  of  change  to  the  momentum  density  of  species  s  is  therefore: 

=  f  d9  93  e~g2/a2  ast{g)  ■  [  dpdce  e29WCe/a2  { g-g')n,  (34) 

Ti2  a6  J  J 
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Let  us  consider  the  last  integral  over  the  scattering  angle.  From  Figure  1,  the  vectors  g,  g' 
in  the  rotated  frame  (£,  rj,  g)  are: 


M 

fcs  \ 

g  =  9  g  =  9  ■ 

0 

;  g'  =  g'  ■  g'  =  g' 

SpSx 

w 

V  c*  / 

Therefore  the  integral  yields: 

J  dQ'(g-g')G(g,Q')  —  2ng  j  dcx5(g,x)g  ~  2irg'  J  dcxcxG(g,x)g 

=  [j-j'{cosx)„,]  g 


(35) 


(36) 


We  must  now  express  the  unit  vector  g  in  the  initial  (x,  y,  z)  frame,  which  is  given  by  (7); 
integration  over  the  ip  variable  leaves  only  one  component,  eg w,  yielding: 


R,  =  -/iw  .  f  dg  g3  e  9~/o?  ast(g)  ■  f  dcece  e29WCe/a2  [g-g'{ cosy)n,]  (37) 

7T2  OiA  J  J 

Using  (20)  to  replace  the  last  integral  and  using  the  normalized  variables  (12)  leads  to: 

Rs  =  --gw  nsntgfe~xJ  dz  z?  e~z  ast(z)  ({1\V\z)  (^z  -  Vz' {cos  x)  (38) 

A  similar  (but  of  opposite  sign)  expression  can  be  obtained  for  the  species  of  type  t,  as  a 
result  of  the  identity  ms(vt— Vt)  =  g(g'  —  g). 

We  can  now  specify  the  type  of  collision.  For  an  £  — y  u  excitation,  we  follow  (15),  to  yield: 
R|  =  --/iw  nsnigfe~x  I  dxx?  eTx  &le(x)  C^1^(v/A^)  (\/x  —  Vx'(cosx)nlJ  (39) 


For  deexcitation,  we  follow  (17): 


Rt  =  -~gwnsnugfe  x 


dx'  (x')a  e  x'  cr\u{x')  (^(VXx')  (\fx’-y/x{cosx)n, 


(40) 


Note  that  the  expressions  (39-40)  are  obtained  in  a  frame  where  w  is  aligned  with  the  z 
direction,  and  corresponds  to  the  change  in  momentum  density  along  that  direction.  Thus, 
it  is  the  component  of  a  force  parallel  to  w,  while  all  components  in  the  transverse  directions 
are  zero,  by  reason  of  symmetry39.  The  components  in  an  arbitrary  rest-frame  must  be 
obtained  by  projecting  w. 
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Since  the  force  density  is  approximately  proportional  to  w,  we  can  group  all  the  other  terms 
into  the  definition  of  a  coefficient,  such  that 


R 1  =  ~Kle(us  -  ut)  Rj  =  +K]e(us  -  ut) 
R t  =  ~kL(us  -  ut)  R|  =  +Kju( us  -  ut) 


where  Kje  and  Kju  are  known  as  resistance  coefficients: 


■77  —A 


Kse  =  3 imsUigpe 

1  2 

Ksu  =  3  gnsnugfe 


-A 


dxx  2e  x  cr^(x)  \/J  —  V/x'(cos  x) 

)  _  r 

da/  (x')5  e-^  <7^(0;)  £^(\/ Xx')  Vx'  —  \/x(cosx)n, 


(41a) 

(41b) 

(42a) 

(42b) 


ft  must  be  pointed  out  that  when  the  collision  is  elastic,  i.e. ,  x  =  a/,  we  recover  the  expression 
of  the  momentum  transfer  cross  section  often  used  in  transport  calculation,  i.e.,  am(x)  = 
<t(1  —  (cosy))  (see  for  example32).  In  the  limit  of  weak  divergence  of  mean  fluid  velocities 
(A  — >  0)  and  isotropic  scattering  (G(x)  —  1  / 47t)  ,  we  have: 


K]t  ~  -jmsnigf  /  dxx2cr\i(x)e~ 


(43) 


Again,  using  the  Klein- Rosseland  relation,  the  excitation  resistance  coefficient  can  be  written 
as: 


*lt  = 


~  1  2  f 00 

Beu(T)  -nnsnegfe~x  / 

J  3  Jo 


dx'x'x^e  x'  C(1)(V A(aa))cr^(rr')  \[x  —  Va/(cosx)n, 


(44) 


As  in  the  case  of  the  zero-th  order  moment,  we  define  the  momentum  exchange  rates  by: 

Kje  =  /msnthile  (45a) 


Kt.  = 


HnsnuKg 


(45b) 


In  the  case  of  weak  divergence  of  mean  fluid  velocities  and  isotropic  scattering,  the  ratio  of 
the  rate  coefficients  for  the  forward  and  backward  processes  is  approximately: 


K 


si 

4. 

^SU 


Beu(T) 


/0°°  dx'e  x'  x'(x*+x')  [l  +  |A(o;*-|-a/)]  cyj;u(a/) 


(46) 


/0°°  dx'e~x'x ,2  [l  +  |  Ax']  a\u(x') 

Note  that  there  is  an  additional  contribution  from  high-order  moment  from  the  expansion. 
This  can  be  seen  by  further  expanding  the  integrand  of  the  numerator: 


0/(0;*+ a/) 


1  +  -A(a;*+a/) 
5 


/2 


=  X 


1  +  -Xx' 


1  *  / 

+  X  X 


1  +  —X(x*  +  2a/) 
5 
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such  that 


Rr  (T\ 

..  jrdx'e  X'x>  i1  +  |A(x*+2a;/)]  ^n(^) 

/0°°  dx'e~x'x'2  [l  +  |Ax']  a\u(x') 

(47) 


Note  that  even  as  A  — >  0,  the  correction  term  does  not  vanish.  Therefore,  there  is  no 
equivalence  between  the  resistance  coefficients  of  the  forward  and  backward  processes  in 
the  limit  A  — >  0.  However,  this  is  perfectly  understandable;  note  that  the  correction  is 
proportional  to  the  energy  threshold,  and  since  kinetic  energy  must  be  removed  from  particle 
s  in  order  to  achieve  excitation,  but  not  for  deexcitation,  there  must  also  be  an  imbalance 
in  the  momentum  exchange  rate.  As  expected,  this  imbalance  vanishes  for  elastic  collisions 
(x*  — y  0)  and  the  rates  are  consistent  with  the  detailed  balance  of  the  mass  exchange.  In  all 
cases,  detailed  balance  is  enforced  through  relation  (26)  at  the  microscopic  level. 


D.  Second-order  moment:  total  energy  density 

The  net  rate  of  change  of  total  energy  of  species  s  can  be  obtained  by  setting  = 
\ms  (v'^  —  v^)  into  equation  (11): 

Qs  =  '  J  d3V*fv*  J  dgg3e~92/a2ast(g)  ■  j  dipdce  e29WCe/a2  (^ms(V2s-v2s))n,  (48) 

Using  the  transformation  defined  in  Appendix  B, 

(v1“vs)  =  h  (g'-g)  '  [V*  +  U— y(g  — w)]  -  (49) 

The  integration  of  the  first  term  in  the  square  bracket  is  zero  since  J  d3V*  V*  fy*  =  0.  One 
can  easily  see  that  the  second  term  in  brackets  is  simply  Rs  •  U  by  comparing  with  (34). 
Similarly,  the  last  term  in  (49)  is  identified  as  —  {writ / M) AeY st.  The  third  term  in  brackets 
involves  the  following  dot  product: 

(g'-g)  '  (g-w)  =  g  ■  g'  -  w  ■  g'  -  g2  +  w  ■  g 

=  gg'  cos  y  —  wg'( g'  •  w)  —  g2  +  wg  cos  6  (50) 


We  can  now  perform  the  averaging  over  the  scattering  angle;  in  particular,  we  have 


j dVL'Q(x)  g'-w  =  (cos x)n,  cos  6 
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so  that 


fdn'g(x)( g'-g)-(g-w)  =  (g' (cos  x}u,-g)<y9-w  cos  9)  (51) 

Thus,  we  obtain,  after  integration  over  (p: 

Qs  =  U  •  Rs  -  ^AeTst  -  fiX~r~e~X  [dgg3ast(g)e-°2/a2  le  (52) 

lvi  7T2a3  J 

where  the  last  angular  integral  is 

Xe  =  \  J dc^9WCe/a2 (cx) o'  9)(g  wco) 

=  W(cx)nl - g )  | g\  J dcee2^a2  -  |  J dcecee2^a2 }  (53) 

=  W(cx)n,-g)g  |c(0)(^Ai)  - 


Therefore,  after  the  change  of  variables  g  — >■  z: 
Qs  — U  •  —  —  AeTst 


+  fI7~Se  X~2  fdZZ‘e  “<T»*W  (v^-V^OOS*),,,)  ((m(V\z)- yC<1|('/V)^) 


(54) 


The  factor  in  front  of  the  integral  can  be  re-arranged  to  yield: 


X‘nsntgfe  A(2 kT) 


Using  the  identity  y(2A;T)  =  jjk(Tt  —  Ts),  the  final  result  has  a  traditional  form: 

Qs  =  u  •  Rs  -  ^A/ffi,  +  Js^k(Tt  -  Ts) 
with  the  thermal  resistance  coefficient  defined  as: 

hi  =  n,ntgte~x  Jdzzie~*aa(z)  (yz -s/z' (cos  x)„)  (<(0,(VV)  -  (lc<1|('tV)') 


(55) 


(56) 


This  result  is  general,  and  we  can  now  make  the  usual  substitutions  for  excitation: 


Ql  =  U  •  Rl  -  +  J]^k(Tt  -  T, 


(57a) 


=  nsnegfe  A  Idxx^e  xo-\e{x)  (^fx-Vx*  (cos  x)S)  fc(0)(^^)  ~  j  (57b) 
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In  the  case  of  deexcitation,  we  can  still  use  the  general  formula  (54),  except  that  in  this  case, 
Ae  =  —  e*.  This  can  be  seen  if  we  start  from  eq.  (48),  which  gives  us: 


Ql  =  -  1*^4"  f  ■  (g'-g)„, 

7 rzar  J 

+  /i7^¥VA  [  d3g'g'al(g')e-^/a2e2^a2(g'~g)  •  (g-w)  (58) 


7T  2  CH° 

+  ^e*T± 
+  M  sv 


Again,  one  can  easily  recognize  the  standard  formulae: 

Qi  =  U  •  Rt  +  ^e*rl  +  Jtjjk{Tt  -  Ts)  (59a) 

Jju  =  nsnugfe~x  Jdx'(x')h-X'alu(x')  (Vx'-^fx( cosx)n,)  ^^(VAx7)  -  y  C(1)(V/Ax7)^ 

(59b) 

which  could  also  be  obtained  directly  from  (55  -  56),  with  the  usual  substitutions  (17). 


We  can  also  express  the  source  term  for  particle  t,  using: 


1  f  12  2\  1  ( (~\t  ms  /\2  /-I  T  ms  \  2 

-m,(vt  -vt)  =  -m,  ((V-— g)  -  (V-— g) 


rns 
2  M 


Mg'  “g  )  -  hV  •  (g'-g) 


(60) 


=  -M g'-g)  •  [V*  +  U  -  7(g-w)]  -  — Ae 


Comparing  with  (49),  we  easily  obtain  (note  the  inversion  of  Ts  and  Tt  in  the  last  term): 

<3?  =  U  ■  Rl  -  +  4  ^k(T,  -  Tt)  (61a) 

Qi  =  U  •  Ri  +  +  4  ^k(T,  -  Tt)  (61b) 

with  R)  =  —  Rs  .  Combining  both  s  and  t  fluids,  the  only  term  remaining  is  the  loss  of 
energy  equal  to  the  energy  gap  between  the  levels,  as  expected.  Note  also  that  this  energy 
loss  is  distributed  to  the  respective  fluids  according  to  the  ratio  of  masses,  such  that  the 
lighter  element  receives  the  major  contribution.  This  is  also  an  expected  result,  similar  to 
the  energy  exchange  due  to  elastic  collisions,  and  due  to  the  kinematics  of  collision. 


In  the  limits  of  near-single  fluid  (A  — y  0)  and  isotropic  scattering,  we  have: 


2A\ 


J It  -  nsn£gf  (  1  -  y  J  J  dxx2  cr]e(x)e 


(62) 
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The  thermal  relaxation  rates  can  be  extracted  similarly: 


(63a) 

(63b) 


At  =  nsndle 

JL  =  nsnujju 


The  ratio  of  the  thermal  relaxation  rates  can  be  written  as: 

;  /0°°  dx'  e~x'  x1  [l  —  |A  +  |A(x*  +  2x')] 


•t 

jU 

jsu 


r  _ 

Blu(T) 

l  +  £* 


cA 

SU 


(64) 


r dx’  e~x‘ x'2 11  -  iA + |w]  <r*„ 

Similar  to  the  case  of  momentum  transfer  rates,  the  correction  terms  do  not  vanish  when 
A  — >•  0  due  to  contribution  from  high-order  moments. 


III.  Numerical  Results 

In  the  following  sections,  we  carry  out  a  numerical  evaluation  and  verification  of  the  ex¬ 
change  rates  derived  in  II B,  II C  and  II D,  for  the  case  of  free  electrons  interacting  with  hy¬ 
drogen  atoms;  these  processes  include  electron- neutral  elastic  collision  and  electron- impact 
excitation  and  deexcitation.  Ionization  and  recombination  are  currently  omitted  and  will 
be  included  in  future  work.  For  comparison  purpose,  we  also  show  the  exchange  rates 
due  Coulomb  collision,  i.e. ,  electron- hydrogen  ion,  which  is  the  dominating  elastic  exchange 
mechanism  for  plasma  with  high  ionization  fraction. 

The  notations  are  slightly  modified  to  better  distinguish  each  type  of  interaction.  We  use 
superscripts  (en),  (ei)  and  (xd)  to  denote  electron-neutral,  electron-ion  (Coulomb),  and 
excitation/dexcitation  (as  a  whole)  collisions,  respectively.  These  grouped  notations  are 
useful,  for  example,  when  looking  at  the  net  momentum  (or  energy)  transfer  due  to  each 
type  of  collision.  The  symbols  'f,  l  are  still  retained  to  indicate  individual  excitation  and 
deexcitation  rates.  For  each  transition  between  two  atomic  states,  we  use  the  convention 
of  indexing  the  final  state  on  the  left,  and  the  initial  state  on  the  right,  i.e.,  For 

example,  is  the  forward  excitation  rate  from  £  to  u,  and  is  the  reverse  process. 

The  energy  levels  and  cross  sections  models  for  atomic  hydrogen  are  given  in  classical  form 
and  summarized  in  Appendix  C 


A.  Reaction  Rates 

All  the  exchange  rates  (mass,  momentum,  energy)  can  be  tabulated  as  a  function  of  two  pa¬ 
rameters:  the  average  thermal  temperature  T,  defined  in  appendix  B,  and  a  non-dimensional 
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parameter  A,  defined  in  eq.  (12),  which  corresponds  to  the  relative  mean  kinetic  energy.  For 
convenience,  we  also  define  an  equivalent  drift  temperature  Tw  =  AT10,  such  that  all  the  rates 
can  be  tabulated  in  terms  of  two  temperatures.  Since  the  mass  ratio  between  the  electron 
and  the  atom  is  very  small  (me  Mh),  we  can  drop  terms  of  order  me/M  and  arrive  at 
the  following  approximations:  /i  ~  me,  T  ~  Te,  Tw  ~  A Te  and  ~gf~ve=  \J There¬ 
fore,  all  the  exchange  rates  for  electron-induced  collisions  (both  elastic  and  inelastic)  can  be 
numerically  evaluated  in  terms  of  the  electron  temperature  Te  and  the  drift  temperature  Tw. 

Figure  2  shows  example  calculations  of  the  zeroth-order  reaction  rates,  defined  in  eqs.  (22), 
(23)  and  (28),  for  electron-impact  excitation  and  deexcitation  between  the  first  three  atomic 
states  of  hydrogen.  These  rates  exhibit  a  similar  trend  for  the  range  of  temperatures  plotted 
here,  that  is,  starting  from  low  temperature,  the  rates  first  increase  reaching  a  plateau  and 
then  decreases  as  temperature  further  increases.  The  value  at  which  the  rate  is  maximum 
is  very  close  to  the  threshold  temperature  of  the  transition.  This  trend  holds  both  in  the 
direction  of  increasing  thermal  Te  (x-axis)  or  drift  temperatures  Tw  (y-axis). 

It  is  clearly  shown  from  Figure  2  that  the  reaction  rates  can  be  significantly  different  from  the 
thermal  limit  when  the  relative  mean  velocity  between  two  fluids  is  significant.  In  particular, 
one  sees  an  increase  of  the  reaction  rate  in  the  low  temperature  regime  where  Te  and  Tw  are 
small  compared  to  the  excitation  temperature  of  the  collision;  this  enhancement  corresponds 
directly  to  the  form  of  the  cross  sections.  Therefore,  one  can  expect  that  significant  deviation 
from  the  thermal  rate  occurs  when  the  mean  kinetic  energy  is  of  the  same  order  as  the 
excitation  temperature.  This  indicates  that  excitation  and  deexcitation  among  low  energy 
states  with  large  threshold  energies  are  more  sensitive  to  the  multifluid  effects.  This  is 
consistent  with  the  prior  statement  we  made  when  examining  the  ratio  of  forward  and 
backward  rates. 

Figure  3  shows  the  forward  and  backward  reaction  rates  of  the  first  transition  between  the 
ground  state  and  the  first  excited  state  as  a  function  of  the  thermal  temperature  for  several 
drift  temperature  values.  Even  at  very  low  thermal  temperature  (Te  ~  0.1  eV),  one  can  have 
significant  excitation  ( w  ~  10~14  m3/s)due  to  a  high  drift  temperature.  In  addition,  Figure  3 
also  shows  that  in  the  limit  A  — >  0,  the  multifluid  rate,  as  formulated  here,  converges  to  the 
expression  for  thermal  limit,  given  in  eq.  (C.3),  as  expected.  Figure  4  shows  the  reaction 
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Figure  2.  Multifluid  reaction  rates  for  electron-impact  excitation/deexcitation  collisions  as  a  func¬ 
tion  of  two  temperatures.  For  better  color  contrast,  all  the  values  smaller  than  minimum  value  of 
the  colormap  are  white  out. 
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Figure  3.  Multifluid  reaction  rates  for  electron-impact  excitation/deexcitation  collisions:  lines  with 
symbols  correspond  to  different  values  of  the  drift  temperature.  The  solid  line  corresponds  to  the 
exact  solution  for  the  thermal  limit  given  by  equation  (C.3). 

rates  as  a  function  of  the  drift  temperature  for  several  thermal  temperature  values;  here 
we  can  identify  a  different  asymptotic  limit  of  the  rate.  In  the  limit  A  — »  oo,  the  electron 
velocity  distribution  function  approaches  the  form  of  a  delta  function  centered  at  the  relative 
mean  velocity  w,  and  the  reaction  rates  approach  the  beam  limit  given  by  eq.  (C.5).  It  must 
be  noted  that  in  the  numerical  integration  of  the  multifluid  rates,  e.g.,  eq.  (22),  the  energy 
grid  x  needs  to  be  refined  near  the  value  of  the  the  mean  kinetic  energy  A  to  avoid  numerical 
error  due  to  the  integration  over  a  delta  function. 

B.  Momentum  and  Energy  Exchange  Rates 

We  now  compute  the  momentum  and  energy  exchange  rates  due  to  both  excitation  and 
deexcitation,  and  compare  with  those  due  to  elastic  collisions  (electron-neutral  and  electron- 
ion).  Recall  from  eq.  (55)  that  the  total  energy  transfer  include  three  terms:  the  first  term 
due  to  work  done  by  friction  in  the  COM  reference  frame,  the  second  term  due  to  thermal 
resistance  (or  thermal  relaxation),  and  the  last  one  due  to  heat  release/absorption  due  to 
chemical  reaction.  There  are  three  different  rate  coefficients  associated  with  each  of  these 
processes,  namely  the  momentum  exchange  rate  k,  thermal  relaxation  rate  j,  and  the  reaction 
rate  w.  In  this  section,  we  will  focus  on  examining  the  momentum  exchange  and  thermal 
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(a)  Excitation  (1  — >•  2) 


(b)  Deexcitation  (1  <—  2) 


Figure  4.  Multifluid  reaction  rates  for  electron-impact  excitation/deexcitation  collisions:  lines  with 
symbols  correspond  to  different  values  of  the  thermal  temperature.  The  solid  line  corresponds  to 
the  exact  solution  for  the  beam  limit  given  by  equation  (C.5). 

relaxation  rates.  The  expressions  derived  in  II C  and  II D  can  be  readily  used  for  the  case 
of  elastic  collisions  by  simply  setting  z*  =  0  and  z'  —  z  —  x.  For  example,  in  the  case  of 
electron-neutral  (en)  collision,  we  have: 

(65a) 
(65b) 

It  can  be  seen  from  the  previous  two  equations  that  in  order  to  compute  the  rate  for  the  case 
of  elastic  collisions,  we  only  need  the  so-called  momentum  transfer  cross  section  a^en'm\x)  = 
u(en)(l  —  (cosy)).  These  cross  sections  are  available  for  a  wide  range  of  neutral  species  due  to 
their  extensive  use  in  calculation  of  transport  properties  (see  for  example13).  The  electron- 
neutral  collision  cross  section  utilized  in  this  work  is  taken  from  Bray  and  Stelbovics34.  For 
inelastic  collisions,  we  need  the  full  DCS,  i.e. ,  both  a(x)  and  G(x,x)  in  ecf  (9)  f°r  each 
process.  These  cross  sections  are  generally  not  available  and  analytical  approximation  is 
needed.  For  Coulomb  collision,  we  use  the  analytical  DCS  from  Rutherford’s  scattering 
formula  with  suitable  cut-off  based  on  the  Debye  length4,  yielding  a  momentum  transfer 
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Figure  5.  Elastic  momentum  transfer  cross  section  for  electrons  in  atomic  hydrogen  computed  with 
the  DCS  from  equation  (68).  The  symbols  are  theoretical  values  from  Bray  and  Stelbovics34. 


cross  section  of  the  form: 
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where  In  A  is  the  well-known  Coulomb  logarithm, 
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A  =  1.24  x  107 

Typically,  In  A  ss  5  —  20;  for  convenience,  we  take  a  constant  value  of  In  A  =  5  for  all  the 
plots  shown  here.  In  this  case,  the  momentum  exchange  and  thermal  relaxation  rates  can 
be  obtained  in  exact  forms28: 
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Due  to  the  lack  of  data  of  the  DCS  for  inelastic  processes,  we  have  used  an  analytical  Born 
scattering  approximation  for  a  Coulomb  screened  potential13,  given  by  the  following  form: 

C  ,  1 
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where  C  is  a  normalization  constant  s.t.  2ir  f\1  Q(e,  x)dcx  =  1.  This  angular-dependent  DCS 
has  been  used  to  compute  momentum  transfer  and  thermal  relaxation  rates  for  both  electron- 
neutral  and  excitation/deexcitation  collisions.  We  note  here  that  the  angular-dependent 
DOS’s  Q  for  excitation  and  deexcitation  have  to  satisfy  detailed  balance30: 

0H€)(£>X)  =  G{i\u)(£’’X)  (69) 

The  above  condition  implies  that  the  probability  for  a  deexcitation  collision  with  an  incident 
energy  e  to  have  a  scattering  angle  of  y  is  the  same  as  that  for  an  excitation  collision  with  an 
incident  energy  £+£*;  therefore,  one  cannot  independently  specify  both  DOS’s  for  the  forward 
and  backward  processes.  Figure  5  shows  a  comparison  of  the  computed  elastic  momentum 
transfer  cross  sections  to  the  result  from  a  direct  close-coupling  calculation  of  Bray  and 
Stelbovics34;  the  agreement  between  the  two  is  excellent.  Using  the  angular-dependent  DCS 
defined  in  eq.  (68),  the  exchange  rates,  e.g.,  n  and  j,  can  be  computed  for  each  bound-bound 
transition  and  summed  over  all  transitions  to  yield  the  total  rates: 

j^{xd)  _ 

Klu\t)nin*  +  4w^)  (70a) 

l  u>£ 

J{  *  =  ^  ^  (j(u\e)nirle  +  i\t\u)nunej  (70b) 

l  U>1 

Based  on  the  total  frictional  and  thermal  resistance  coefficients,  we  can  extract  average 
momentum  transfer  and  thermal  relaxation  rates  for  all  excitation/deexcitation  processes  as 
follows: 


K<,xd)  —  menennK^x<^  (71a) 

jOrf)  =  nennj{xd)  (71b) 

where  rin  =  l^ie  total  atomic  number  density  (summation  over  levels).  Note  that 

according  to  our  definitions,  the  exchange  rates  K('xd'>  and  contain  terms  designating  the 
population  of  the  excited  states,  e.g.,  rik/nn.  For  comparison  purpose,  these  average  rates  are 
calculated  by  assuming  a  Boltzmann  distribution  of  the  atomic  states,  i.e. ,  riik  =  nn9ke  ^  — — 
and  Zn  =  9ke ~Ek/kTB.  Once  can  see  that  in  this  case  the  population  of  the  excited 

states  is  effectively  replaced  by  a  Boltzmann  distribution  characterized  by  a  temperature  Tb- 
This  step  is  only  done  for  the  comparison  shown  below.  In  a  CR  calculation,  the  detailed 


23 


population  of  the  atomic  states  (equilibrium  or  not)  is  known  and  the  exchange  rates  are 
computed  for  each  transition  as  specified  in  eq.  (70). 

Figure  6  shows  a  comparison  between  the  momentum  exchange  and  thermal  relaxation 
rates  for  three  different  processes:  electron-neutral,  electron-ion  (Coulomb),  and  excita¬ 
tion/deexcitation  collisions.  For  clarity,  the  rates  due  to  Coulomb  collision  are  only  shown 
for  value  of  Te  >  1  eV.  Two  important  observations  can  be  deduced  from  this  plot.  Firstly, 
when  the  atoms  are  cold,  i.e. ,  Tg  is  low,  the  inelastic  exchange  rates  are  much  smaller  than 
elastic  ones.  This  is  due  to  the  fact  that  when  Tg  is  low,  only  collisions  between  the  ground 
state  and  a  first  few  excited  states  are  significant;  the  rates  for  these  transitions  are  low  com¬ 
pared  to  the  others  due  to  larger  energy  threshold.  As  the  atoms  are  being  excited  and  Tg 
increases,  transitions  among  highly  excited  states  become  significant,  leading  to  an  overall 
increase  in  the  total  rates.  These  rates  eventually  exceed  those  due  to  elastic  collisions  as 
can  be  noticed  in  region  (I)  in  Figure  6  for  the  dash-dotted  and  dotted  lines.  Secondly,  the 
rates  due  to  inelastic  processes  tends  to  have  a  slower  drop-off  at  high  temperature  compared 
to  the  elastic  rates,  which  suggests  that  at  sufficiently  high  temperature,  the  main  momen¬ 
tum  and  energy  transfer  mechanisms  (region  (If)  in  Figure  6)  are  due  to  inelastic  collisions. 
Thus,  we  are  led  to  the  important  conclusion  that  one  cannot  neglect  momentum  and  energy 
transfer  due  to  inelastic  collisions.  The  only  justification  for  neglecting  these  terms  is  when 
the  atoms  are  cold  and  thermal  temperature  is  low,  but  both  of  these  conditions  will  not  be 
realized  in  most  practical  systems. 

C.  Verification 

The  accuracy  of  the  derived  formulas  of  the  exchange  source  terms  are  verified  against  direct 
evaluation  of  the  full  transfer  integral  (3)  over  six  dimensional  space  using  Monte  Carlo 
method.  The  procedure  for  the  Monte  Carlo  integration  is  as  follows:  (1)  sample  different 
pair  of  particles  (one  atom  and  one  electron)  from  two  different  Maxwellian  distributions, 
(2)  compute  the  exchange  rate  due  to  each  sample  pair,  and  (3)  accumulate  these  rates.  The 
sum  of  these  rates  will  follow  the  correct  probability  distribution  function  of  the  samples. 
For  brevity,  we  only  show  an  example  calculation  of  the  zeroth-order  reaction  rate  for  an 
excitation  and  deexcitation  from  levels  2  and  5.  Figure  7  shows  a  detailed  comparison 
between  the  rate  expressions  obtained  from  eqs.  (22)  and  (23)  against  Monte  Carlo  results 
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(a)  Momentum  exchange 


(b)  Thermal  relaxation 


Figure  6.  Comparison  of  the  momentum  exchange  rates  n  and  thermal  relaxation  rates  j  due  to 
elastic  and  inelastic  collisions  (excitation  and  deexcitation).  Different  colors  of  the  line  indicate 
different  values  of  the  drift  temperatures  in  eV:  Tw  =  0.1  (blue),  1  (green)  and  10  (red)  eV.  Solid 
lines  indicate  rates  due  to  e_-H  elastic  collisions,  solid  lines  with  circles  indicate  e--H+  (Coulomb) 
collisions,  and  three  sets  of  broken  lines  are  used  for  the  inelastic  collisions.  The  rates  for  inelastic 
collisions  are  defined  in  eq.  (71)  and  are:  dashed,  dash-dotted  and  dotted  lines  for  the  case  of 
Tb  =  0.5, 0.7  and  0.8  eV,  respectively. 

with  an  excellent  agreement.  Similar  agreement  is  obtained  with  other  sampled  transitions, 
giving  us  complete  confidence  in  the  accuracy  of  our  derivation  of  the  multifluid  rates  from 
kinetic  theory. 

D.  Zero-dimensional  Calculations 

We  conducted  zero-dimensional  (0D)  calculations  for  a  constant-volume  (isochoric)  system 
to  support  our  findings  in  sections  III  A  and  III  B,  which  include  the  following:  (1)  reaction 
rates  can  be  enhanced  when  the  relative  mean  drift  velocity  is  significant,  and  (2)  momentum 
transfer  and  thermal  relaxation  due  to  inelastic  collisions  are  non- negligible.  The  system 
contains  hydrogen  atoms  and  a  small  fraction  of  free  electrons41.  The  governing  equations 
are  described  in  appendix  D,  and  the  resultant  system  of  ordinary  differential  equations  is 
solved  using  the  Radau5  method  of  Hairer  and  Wanner36,  which  is  ideally  suited  for  stiff 
problems. 
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(a)  Excitation  (2  — »  5)  (b)  Deexcitation  (2  •<—  5) 


Figure  7.  Comparison  of  zeroth-order  reaction  rates  with  Monte  Carlo  integration  of  the  full  transfer 
integral.  Lines  are  from  our  derived  expressions.  Symbols  are  Monte  Carlo  results.  For  each  line, 
the  ratio  of  the  drift  and  thermal  temperature  A  is  fixed  and  the  values  are  shown  in  the  figure. 


number  density 

temperature 

atomic 

states 

rii  =  0.9 nt  for  i  =  1 

rii  =  10~15nt  otherwise 

0.3  eV 

electron 

ne  =  0.1  nt 

2  eV 

Table  I.  Initial  conditions  of  OD  test  cases.  For  all  cases,  the  total  atomic  density  n*  is  1020  nr  3. 

The  initial  conditions  for  the  number  densities  and  temperatures  of  the  atoms  and  free 
electrons  are  summarized  in  Table  I.  Initially,  all  the  atoms  are  at  the  ground  state  (denoted 
by  i  —  1  in  table  I)  with  a  translational  temperature  of  0.3  eV.  The  atoms  are  assumed 
to  be  at  rest,  i.e. ,  their  mean  velocity  is  zero.  The  free  electrons  have  a  temperature  of  2 
eV,  and  their  mean  velocity  is  varied  to  demonstrate  the  multifluid  effects.  In  all  the  test 
cases,  we  include  the  first  10  atomic  levels  of  hydrogen  according  to  the  model  described  in 
Appendix  C. 

In  the  first  test,  we  perform  the  calculation  with  various  initial  mean  velocities  of  the  elec¬ 
trons,  or  equivalently,  the  drift  temperatures  Tw.  Figure  8  shows  the  time  evolution  of  the 
number  densities  of  the  excited  states  during  the  isochoric  heating  process  for  two  cases: 
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Figure  8.  Number  density  of  excited  states  during  a  zero-dimensional  chemistry  test.  Solid  lines 
correspond  to  the  solutions  where  the  relative  mean  velocity  is  very  small  (Tw  =  0.01  eV).  Dashed 
lines  correspond  to  the  solutions  with  a  large  relative  mean  velocity  (Tw  =  10  eV). 

Tw  =  0.01  and  10  eV;  the  former  corresponds  to  a  single  fluid  calculation  of  a  thermal  bath 
with  a  warm  electron  population,  and  the  latter  corresponds  to  a  situation  where  an  electron 
beam  is  injected  to  the  system.  One  can  clearly  see  from  Figure  8  that  there  is  an  enhance¬ 
ment  to  the  excitation  process,  indicated  by  an  early  increase  in  the  population  of  excited 
states,  due  to  the  presence  of  a  non-zero  mean  velocity  of  the  electrons.  The  same  argument 
can  be  made  from  Figure  9,  which  shows  the  time  evolution  of  the  temperatures  for  three 
cases  of  different  initial  Tw.  In  this  plot,  the  Boltzmann  temperature  Tb  indicates  the  degree 
of  excitation  of  the  atom.  It  must  be  pointed  out  that  the  enhancement  in  excitation,  how¬ 
ever,  persists  on  the  time  scale  of  the  momentum  relaxation  process,  which  is  indicated  by 
a  drop  in  Tw  at  approximately  3  x  10-7  sec  as  shown  in  the  bottom  plot  of  Figure  9.  After 
this  time,  the  momentum  of  the  electrons  is  completely  absorbed  by  the  atom,  signifying  a 
change  to  single  fluid  kinetics.  In  all  the  test  cases,  the  excitation  proceeds  at  time  scale 
much  smaller  than  the  resolution  of  the  figures,  i.e.,  Tb  approximately  goes  from  0  to  0.7 
eV  in  10~9  sec.  As  mentioned  before,  when  Tb  is  sufficiently  large,  the  momentum  exchange 
and  thermal  relaxation  rates  from  inelastic  collisions  cannot  be  neglected. 

In  the  second  test,  we  specifically  identify  the  effect  of  inelastic  collisions.  Figure  10  shows 
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Figure  9.  Time  evolution  of  the  temperatures  for  several  test  cases  with  different  initial  drift 
velocities:  solid  line  ( Tw  =  0.01  eV),  dotted  line  ( Tw  =  1  eV),  and  dashed  line  ( Tw  =  10  eV).  The 
bottom  plot  shows  the  evolution  of  the  drift  temperature  Tw. 

a  comparison  of  the  temperature  evolution  for  two  test  cases,  both  with  a  very  small  initial 
drift  temperature  Tw  =  0.01:  the  solid  lines  correspond  the  solution  with  both  elastic  and 
inelastic  exchanges,  and  the  dashed  lines  to  the  solution  without  inelastic  exchanges.  In 
this  case,  the  friction  is  negligible  since  Tw  &  0,  and  there  are  only  thermal  relaxation 
and  heat  release  from  reaction.  It  is  clearly  shown  that  inelastic  collisions  do  contribute  to 
the  total  energy  transfer  between  the  electrons  and  atoms,  leading  to  a  faster  temperature 
equilibration  Tn  —  Te  —  Tb-  It  is  interesting  to  point  out  here  that  while  Ts  —  Tb  equilibration 
is  due  to  the  heat  release/absorption  term,  Te  —  Tn  equilibration  is  only  due  to  the  thermal 
relaxation  term.  Figure  10  indicates  that  there  is  a  faster  Te  —  Tn  equilibration  when  inelastic 
thermal  relaxation  is  included.  This  result  is  quite  intriguing,  since  this  effect  is  normally 
neglected  in  most  of  single  fluid  multi-temperature  calculations. 

Figure  11  shows  the  results  for  the  same  two  test  cases  but  now  with  an  initial  drift  temper¬ 
ature  Tw  =  10  eV.  One  can  make  the  same  argument  that  inelastic  collisions  further  enhance 
the  temperature  equilibration  process.  However,  it  is  worthwhile  to  point  out  that  during 
the  time  period  10~8  <  t  <  10~',  the  electrons,  although  being  decelerated  due  to  friction, 
are  also  heated  due  to  the  work  done  by  the  same  force  (first  term  on  the  right  hand  side  of 
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Figure  10.  Time  evolution  of  the  temperatures  for  the  case  with  (solid)  and  without  (dashed)  the 
exchange  terms  (momentum  and  energy)  due  to  inelastic  collisions.  Initially,  Tw  =  0.01  eV. 

eq.  (D.5a)).  This  thermal  heating  term  contains  contribution  from  both  elastic  and  inelas¬ 
tic  collisions,  which  explains  why  the  electron  temperature  drops  faster  towards  equilibrium 
when  inelastic  exchanges  are  neglected.  Finally  we  should  emphasize  that  the  thermal  equi¬ 
librium  between  all  components  (Th,  Te,  Tv.)  can  be  obtained  purely  from  inelastic  collisions, 
as  a  result  of  properly  accounting  for  detailed  balance  in  our  model;  elastic  collisions  are 
not  required  to  achieve  equilibrium,  but  of  course  are  needed  to  obtain  the  correct  rate  of 
relaxation. 

IV.  Concluding  Remarks 

We  have  presented  a  model  for  inelastic  collisions  for  electronic  excitation  and  deexcitation 
within  the  context  of  a  multifluid  description  of  a  plasma.  The  model  is  rigorously  derived 
from  kinetic  theory  and  is  applicable  to  any  multifluid  plasma,  irrespective  of  the  mass 
ratio,  and  strictly  obeys  the  detailed  balance  principle.  The  appropriate  mass  transfer  rates 
and  momentum  and  thermal  resistance  coefficients  are  derived,  and  are  found  to  satisfy  the 
proper  asymptotic  limits,  such  that  in  the  limit  of  vanishing  energy  gap,  the  well-known 
expressions  for  elastic  collisions  are  recovered. 

Numerical  evaluations  of  the  multifluid  rates  are  carried  out  for  a  two-fluid  electron- hydrogen 
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Figure  11.  Time  evolution  of  the  temperatures  for  the  case  with  (solid)  and  without  (dashed)  the 
exchange  terms  (momentum  and  energy)  due  to  inelastic  collisions.  Initially,  Tw  =  10  eV.  The 
bottom  plot  shows  the  evolution  of  Tw. 


plasma  using  Bohr  model  for  the  energy  levels  and  semi-classical  cross  sections.  Several  nu¬ 
merical  tests  were  performed  in  a  virtual  (zero-dimensional)  test  cell,  and  both  the  known 
thermal  and  beam  limits  were  correctly  recovered  from  the  model.  We  also  found  that  in 
some  plasma  conditions  of  interest,  the  contribution  of  the  inelastic  collisions  to  the  resis¬ 
tance  coefficients  is  significant,  contrary  to  the  usual  assumptions  made  in  current  multifluid 
models.  While  this  work  is  focused  on  two-body  reactions,  extension  to  ionization  and  re¬ 
combination  processes  is  currently  under  examination,  and  will  be  presented  in  a  future 
article. 


Acknowledgments 

The  authors  thank  Drs.  Ann  Karagozian,  Justin  Koo,  Robert  Martin,  and  Owen  Smith 
for  stimulating  discussions.  Research  was  supported  by  the  Air  Force  Office  of  Scientific 
Research  (AFOSR),  grant  numbers  12RZ06COR  (PM:  Dr.  F.  Fahroo)  and  14RQ05COR 
(PM:  Dr.  J.  Marshall). 


30 


Appendix  A  Collision  kinematics 


Let  us  consider  an  inelastic  collision  between  two  particles  s  and  t,  such  that  the  post-collision 
particles  can  have  modified  internal  states.  The  process  is  formally  described  as  the  relation 

s(vs)  +  t(vt)  -»•  s'(v'a)  +  (A.l) 


Note  that  only  two  particles  are  produced  by  the  collision.  The  initial  velocities  are  vs,vt. 
The  mean  fluid  velocity  is  u,  such  that  u  =  (v)  =  f  d3vvf(v)  and  a  thermal  velocity 
c  =  v  —  u.  By  definition,  we  also  have  (c)  =  0. 

The  collision  can  be  transformed  to  the  center  of  mass  (COM)  reference  frame,  moving  with 
velocity  V  with  respect  to  the  LAB  frame.  Similarly,  we  can  also  define  a  mean  velocity 
of  this  COM  frame  as  U.  The  subsequent  Galilean  transformations  yield  the  following 
definitions: 


V  = 

msys  +  mtyt 

g  =  vs-vt 

M 

u  = 

msus  +  rntut 

w  =  us  -  ut 

M 

where  M  =  ms+mt.  The  inverse  transformation  yields: 


vs  =  V 


mt 

M 


g 


,,  ma 

vf  =  V - g 

*  MB 


u,  =  U 


mt, 

M 


-w 


TT  ms 

uf  =  U - — w 

M 


(A. 2a) 
(A. 2b) 


(A. 3a) 
(A. 3b) 


Mass  conservation  imposes  the  relation  ms-\-mt  =  M  =  m's-\-m't.  For  the  case  of  two-body 
processes  such  as  excitation  of  internal  states,  the  masses  are  individually  conserved,  i.e. 
m's  =  ms,m't  =  mt.  Expressed  in  the  COM  frame,  momentum  and  energy  conservation  yield, 
respectively: 


MV  —MY'  (A. 4a) 

l-MW2  +  i«52  =i  MV'2  +  +  Ae  (A. 4b) 

where  n  —  msmt/M.  Therefore,  we  have  the  following  constraints: 

o  A  f 

V  =  V'  and  g2  =  g '2  H -  (A. 5) 
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For  an  excitation  between  two  atomic  levels,  the  transferred  energy  is  a  fixed  value  Ae  = 
e*  >  0,  the  energy  gap  between  levels.  For  a  deexcitation,  we  use  the  same  formulation  as 
above  (i.e.  post-collision  variables  indicated  by  a  prime),  but  in  this  case,  Ae  =  —  e*  <  0.  In 
the  limit  Ae  — >•  0,  the  collision  is  elastic. 


Appendix  B  Separation  of  variables 

Consider  the  Maxwellian  velocity  distribution  functions  (VDF)  of  each  particle  type,  nor¬ 
malized  to  unity,  e.g.  (recall  that  c  =  v  — u): 


f- 


fs(vs)  =  (  0  ]  exp 


V  2irkTs 


msc2s 

2kTq 


(B.l) 


and  similarly  for  ft.  The  averaging  over  initial  states  will  yield  a  product  of  these  two 
distributions: 

3  3 

(  m  V 


/a(va)/t(vt)  = 


^  2nkTs )  \  2nkTt ) 

where  the  argument  of  the  exponential  function  is,  from  (A. 3): 

2 


exp  [A] 


A  = 


ms 
2  kTq 


V-U  +  — (g- w) 


+ 


mt 

2kTt 


777  1  2 

V  —  U - -(g  —  w) 


(B.2) 


(B.3) 


Following  Burgers28,  this  expression  can  be  simplified  with  an  appropriate  transformation  of 
variables;  since  the  basic  procedure  will  be  used  elsewhere,  we  describe  it  below.  First,  we 
define  the  following  variables 


Pp  = 


mp 

2kT,n 


g  =  g  -  w 


(B.4) 


such  that 


A  =  ps  (V-U)  +  ^g]  +  fit  [(V— U)  —  -jjg 


=  (&+/3t)(V-U)2  + 

Dehne  now 

and  comparing  the  expression 


lt  n  m2s 

'’NP+l3tNP 


+  g2  +  2 


v*  =  V-U  +  7g 


(B.5) 


(B.6) 


(ps+pt)v*2  =  ( ps+pt)(\-uy  +  (ps+Pt)r  g2  +  27G9.+/3t)(v-u)  •  g 


(B.7) 
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with  (B.5),  we  can  choose  the  appropriate  value  of  the  coefficient  7  to  eliminate  the  dot 
product  from  A: 


1  ( pP+_ 

7  A+A  Vs  M  PtM 

We  then  obtain  complete  separation  of  variables: 


A=(ps+pt)V*2  + 


mt  „  mz 


PsT7o+Pt 


—  -  ft  —V 


'  M2  M2  Ps+Pt  \  M 
The  term  in  brackets  is  easily  simplified: 

r  1  =  PaPt 

Ps  +  Pt 

We  can  now  define  effective,  average  temperatures: 

o  ,  a  _  rns  rnt  M  rnsTt  +  mtTs  _ 

Ps  +  Pt  —  —  + 


M) 


(B.8) 


(B.9) 


(B.10) 


2  kTs  2  kTt  2k  MTsTt 

PsPt  P  M  TsTt  _  p 


M 

2  kT* 


and  7  becomes: 


Ps  +  Pt  2 A;  TsTt  msTt+mtTs  2kT 


P  Tt-Ts  Tt  -  Ts 

7  =  —  — - —  =  P- 


M  T  r"  msTt  +  mtTs 
To  summarize,  we  have  performed  the  following  change  of  variables: 


v*  =  V  -  U  +  p 


Tt-Ta 

msTt+mtTs 


T*  =  M- 


TsTt 


g  =  g  -  w 

msTt+mtTs 


T  = 


(B.lla) 

(B.llb) 

(B.12) 

(B.13a) 

(B.13b) 


msTt+mtTs  M 

These  are  the  same  expressions  found  in28  (pp.  45-46)  (with  an  occasional  change  of  naming 
convention)  for  which  it  is  easy  to  verify  that  the  Jacobian  of  the  transformations  is  unity, 
i.e. 

d3vsd3vt  =  d3Vd3  g  =  d3V*d3  g  (B.14) 

Furthermore,  we  note  that: 

M  Y  f  p  5  ^ 


™s\*  (  mt  \  *  _  (O  o  N  _  (O  ,  a  A  (  PsPt  X  " 

2kTq  )  \  2 kTt  I  ~  ^  ^ 


Ps+Pt 


The  product  of  two  distributions  can  now  be  written  as: 


2  kT*  J  \2kT 


(B.15) 


r,  ■  /<  =  (  M 


2-jikT* 


exp 


MV 


*2 


2  kT* 


p 


2nkT 


exp 


FT 
2  kT 


=  /*(V*)-/(  g)  (B.16) 
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All  subsequent  expressions  can  now  be  simplified  with  this  separation  of  variables.  For 
example,  any  operator  G  that  depends  only  on  variables  expressed  in  the  COM  frame  (g,  g'), 
we  have: 

[  d3vsd3vtfsftG( g,g')  =  [  d3V*r  (V*)-  [  d3gf(g)G(g,g')  (B.17) 


=  1 


Note  that  this  procedure  applies  equally  well  for  elastic,  excitation  and  deexcitation  collisions, 
and  that  no  approximations  have  been  made  on  the  mass  ratio.  Furthermore,  since  the 
averaging  over  initial  states  only  involves  the  distribution  functions  for  the  s  and  t  particles, 
we  have  not  necessarily  assumed  that  the  final  products  s'  and  t'  belong  to  the  same  fluid 
as  the  initial  particles. 

Appendix  C  Atomic  data  and  cross  section  models 

The  atomic  states  of  the  Hydrogen  atom  are  listed  as  a  function  of  their  principal  quantum 
number  {n)  only,  following  the  Bohr  atomic  model;  the  splitting  of  states  with  respect  to 
orbital  and  spin  numbers  is  ignored,  and  all  states  have  a  degeneracy  gn  =  2 n2.  The  states 
number  from  n  =  1  to  oo  and  we  consider  a  finite  number  of  states  n  —  1, . . . ,  M  <  oo  before 
reaching  the  ionization  limit42.  In  this  simplified  model,  the  energy  of  each  state  is  given 
as  En  =  IH  (1  —  1/n2),  as  measured  from  the  ground  state  (Ei  =  0),  and  we  will  denote  by 
In  =  IH  (1/n2  —  1/M2)  ~IH/n2  the  energy  required  for  ionization  of  level  n. 

The  classical  form  of  the  cross  section  for  energy  exchange  between  a  free  electron  and  the 
atom  (Hydrogen)  is  used5.  For  an  excitation  collision  from  level  £  to  level  u  >  £,  the  cross 
section  takes  the  form: 


a 


t 

(u\£) 


[x)  =  (47ra2)(3f&) 


lh 


( X  -  xlu) 


(C.l) 


kTej  x£ux^ 

where  a0  is  the  Bohr  radius,  x  is  the  nondimensional  incident  energy  of  the  electron,  x£u  = 
(Eu  —  Ei)/kTe  is  the  energy  gap  between  £  and  u,  and  is  the  oscillator  strength: 

32  1  1  1 


Uii  — 
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In  the  thermal  (singlefluid)  limit  (A  — >  0),  the  reaction  rate  can  be  obtained  in  an  exact 
form: 


^Uf)  -  (4vra2)(3ffa)ue 


m 


(C.3) 
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where 


J  Iptu 


El  (Xiu)  and  Ei(a;)  = 


r°°e -y 


(C.4) 


Here,  ve  is  the  mean  thermal  electron  velocity  and  Ei  is  the  exponential  integral.  On  the 
other  hand,  in  the  beam  limit  (Te  — y  0),  the  reaction  rate  takes  the  form: 


where 


~  A)  =  (4vra0)(3fte)n 


y^r  (A  —  xlu) 
2  xiu\2 


■  (Iz. 
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Appendix  D  Rate  equations 

We  describe  here  the  governing  equations  for  the  zero-dimensional  simulation,  which  describe 
the  time  evolution  of  a  constant-volume  system  during  a  thermochemical  relaxation  process. 
The  following  system  of  rate  equations  is  considered: 


[mPnPuP)  —  —  Up 


un)(/i(en)  +  K(xd)) 


—  (mnnnun)  =  +(ue  -  u n){K(en)  +  K(xd)) 


=  —  U  •  Up  —  u, 


t)(/i(en)  +  I&d)) 


-^k(Te-Tn)(J^  +  jW)- 
mn 

=  +U-(ue-U  n)(K^  +  K^) 

+  —k{Te  -  Tn)(J{en)  +  J(xd) ) 


(D.la) 

(D.lb) 

(D.lc) 

(D.ld) 


(D.le) 


(D.lf) 


where  nn  =  ^2ien  ni  and  -Ee(n)  is  the  total  energy  of  the  fluid: 

-S'e(n)  ^ne(n)kTe(n)  H“~^'e(n)^'e(n)^-e(n)  *  ^-e(n) 


(D.2) 
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where  £etn)  is  defined  as  the  thermal  energy  of  the  fluid.  The  exchange  source  terms  can  be 
decomposed  into  two  parts:  elastic  and  inelastic  collisions  between  the  electrons  and  atoms. 
The  momentum  and  energy  exchange  rates  due  to  elastic  and  inelastic  collisions  are  defined 
in  eqs.  (65)  and  (70).  The  remaining  terms  are: 


-p(xd) 

1  i 


^(xd) 


[  nineW\u\i)  +  nunett\i\u) 

u>i 


77/2 77/ e TliTleZU 


EE  ninewju^e£u  —  nuneu 

i  u>£ 


(D.3) 

(D.4) 


where  e\u  =  Eu  —  E(_.  Comparing  eqs.  (D.lc)-(D.lf)  and  (D.2),  we  can  also  write  conservation 
equations  for  the  thermal  energies  of  electrons  and  atoms: 

rip 

_J1  =  -me(u  -  Ue)  •  (ue  -  u n)(K^  +  K M) 

-  —k(Te  -  Tn)(J{en)  +  J(xd) )  -  E^xd)  (D.5a) 

mn 

dr 

=  +me(u  -  un)  •  (ue  -  u n)(K^  +  KW) 

+  —k{Te  -  Tn)(J (en)  +  J{xd) )  (D.5b) 
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